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Abstract 

Surfaces serve as highly efficient catalysts for a vast variety of chemical reactions. Typically, 
such surface reactions involve billions of molecules which diffuse and react over macroscopic areas. 
Therefore, stochastic fluctuations are negligible and the reaction rates can be evaluated using rate 
equations, which are based on the mean-field approximation. However, in case that the surface 
is partitioned into a large number of disconnected microscopic domains, the number of reactants 
in each domain becomes small and it strongly fluctuates. This is, in fact, the situation in the 
interstellar medium, where some crucial reactions take place on the surfaces of microscopic dust 
grains. In this case rate equations fail and the simulation of surface reactions requires stochastic 
methods such as the master equation. However, in the case of complex reaction networks, the 
master equation becomes infeasible because the number of equations proliferates exponentially. 
To solve this problem, we introduce a stochastic method based on moment equations. In this 
method the number of equations is dramatically reduced to just one equation for each reactive 
species and one equation for each reaction. Moreover, the equations can be easily constructed 
using a diagrammatic approach. We demonstrate the method for a set of astrophysically relevant 
networks of increasing complexity. It is expected to be applicable in many other contexts in which 
problems that exhibit analogous structure appear, such as surface catalysis in nanoscale systems, 
aerosol chemistry in stratospheric clouds and genetic networks in cells. 

PACS numbers: 05.10.-a,82.65.+r,98.58.-w 
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I. INTRODUCTION 



The catalysis of chemical reactions by surfaces plays a crucial role in a vast range of 
physical, chemical and biological systems. In many cases, the surfaces are of macroscopic 
dimensions and the reactants appear in large quantities. Under these conditions, the law of 
large numbers applies and fluctuations in the surface concentrations of reactants and their 
reaction rates become negligible. As a result, these reactions can be analyzed using rate 
equation models that incorporate the mean-field approximation and ignore fluctuations. 

Consider the case of a surface, which is partitioned into microscopic domains, that are 
disconnected, namely reactants cannot diffuse between them. The population of reactive 
atoms and molecules in each domain becomes small and fluctuations become important. As 
a result, rate equations fail and the simulation of these reactions requires stochastic methods 
such as direct integration of the master equation 1, or Monte Carlo (MC) simulations 

. While MC simulations require the accumulation of statistical data over long times, 
the master equation provides the probability distribution from which the reaction rates can 
be obtained directly. In certain cases, the master equation can be solved using a generating 
function , llO|, • The set of coupled ordinary differential equations is then trans- 



formed into a single partial differential equation for the generating function. This equation 



can solved numerically and in a few cases can also be solved analytically [12j, [13] . The master 
equation can also be approximated using the Fokker-Planck equation This is a partial 
differential equation in which the population sizes of the reactive species are represented 
by continuous variables 15) . The master equation and related methods described above 
are useful for simple reaction networks, which involve few reactive species. However, as the 
number of reactive species increases, the number of equations in the master equation quickly 
proliferates. This makes the master equation infeasible for complex networks including a 
large number of reactive species. 

The recently introduced multiplane method provides a dramatic reduction in the number 
of equations Iff]. In this method, the reaction networks are described by graphs, where each 
node represents a species and each edge represents a reaction. Typically, these networks are 
sparse, namely most species react only with few other species. In the multiplane method one 
breaks the network into a set of fully connected subnetworks (cliques). Lower- dimensional 
master equations are constructed for the marginal probability distributions associated with 
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the cliques, with suitable couplings between them. This enables the simulation of large 
networks much beyond the feasibility limit of the master equation. However, it turns out 
that the construction of the multiplane equations for complex networks is difficult. 

In this paper we present a method for stochastic simulations of surface reaction networks, 
which is based on moment equations. This method exhibits crucial advantages over the 
multiplane method 17|. The number of equations is further reduced to one equation for 



each reactive species (node) and one equation for each reaction (edge). Thus, for typical 
sparse networks the complexity of the stochastic simulation becomes comparable to that of 
the rate equations. Unlike the master equation (and the multiplane method) there is no need 
to adjust the cutoffs, namely the same set of equations is used under all physical conditions. 
Moreover, for any given network the moment equations can be easily constructed using a 
diagrammatic approach. This enables to automate the construction of the set of equations 
- a feature which is essential in the case of complex networks. Furthermore, the moment 
equations are linear in terms of the moments. Thus, the stability and convergence properties 
can be easily controlled and the steady state solution can be obtained by standard algebraic 
methods. 

Below we consider some examples of surface reaction systems in which fluctuations play 
an important role. In surface catalysis systems, the surface is often partitioned into facets 
on nanometric dimensions, with little diffusion between the facets. Under these conditions 
the number of reactive atoms and molecules residing on a facet is small and fluctuations are 



strong 



. 3,Q,y, 
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23| . Another important example of chemical reaction networks 
that require stochastic analysis appears in the field of interstellar chemistry. Some chemical 
reactions in interstellar clouds take place on the surfaces of dust grains 
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1261,1271. 



as well as reaction networks 



es. Unlike gas phase reactions in cold 



These include molecular hydrogen formation 
that form ice mantles and certain organic molecu 
clouds that mainly produce unsaturated molecules 
hydrogen-addition reactions that result in saturated, hydrogen-rich molecules, such as water 
HqQ), ammonia (NH3), methane (CH4), formaldehyde (H2CO) and methanol (CH3OH) 



32j, surface processes are dominated by 



33. 



34J. In particular, recent experiments have shown that methanol cannot be efficiently 



produced by gas phase reactions 35|]. On the other hand, there are indications that it can 



be efficiently produced on ice-coated grains 



Therefore, the ability to perform efficient 



simulations of chemical reactions on interstellar grains is of great importance. 
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Due to the submicron size of the grains and the low gas density, the populations of 



reactive species 



suitable 
equation 
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)er grain are small and strongly fluctuate. Therefore, rate equations are not 
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40l | and stochastic methods such as the direct integration of the master 



4l| or MC simulations 42J are required. As discussed above, these methods 
apply in the case of small networks but become infeasible for large networks 43j, l44j ]. 

Here we demonstrate the moment-equation method for a set of astrophysically relevant 
networks of increasing complexity, culminating in the network that describes the production 
of methanol. It is shown that in spite of the small number of equations and the ease of their 
construction, the accuracy of the results is not compromised. The results of the moment 
equations are in excellent agreement with the master equation and coincide with the rate 
equations in the limit of large grains. 

The paper is organized as follows. In Sec. II we consider the production rate of molecular 
hydrogen using the rate equation, master equation and moment equation methods. In Sec. 
Ill we consider the simulation of complex reaction networks using the moment equations. 
The construction of the set of moment equations using diagrammatic methods is presented 
in Sec. IV, followed by a summary in Sec. V. 



II. MOLECULAR HYDROGEN FORMATION 



A. Rate Equations 

Consider the reaction H+H — >H 2 that takes place on the surfaces of interstellar dust 
grains. Hydrogen atoms from the gas phase collide and stick to the surface of a dust grain. 
They diffuse between adsorption sites on the surface and when two of them encounter each 
other they recombine and form an H2 molecule, which desorbs into the gas phase. 

For simplicity, we assume that the grains are spherical and denote their diameter by 
d. The cross-section of a grain is a = nd 2 /4 and its surface area is ird 2 . The density of 
adsorption sites on the surface of a grain is denoted by s (sites cm -2 ). Thus, the number of 
adsorption sites on the grain is S = nd 2 s. Each grain is exposed to a flux F = nva (s -1 ) of 
H atoms, where n (cm -3 ) is the density of H atoms in the gas phase and v (cm s -1 ) is their 
average velocity, which is determined by the gas temperature. It is also useful to define the 
flux in units of monolayers (ML) per second, namely / = F/S (ML s _1 ). Clearly, this flux 
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is independent of the grain size. 

The desorption rate of H atoms from the grain is given by W = z/-exp(— Ei/ksT), where 
v is the attempt rate (standardly taken to be 10 12 s -1 ), E\ is the activation energy barrier 
for desorption of an H atom, and T (K) is the surface temperature. The hopping rate of 
adsorbed atoms between adjacent sites on the surface is a = v • exp(— E^/fceT), where Eq is 
the activation energy barrier for hopping of the atoms. 

Molecular hydrogen formation on astrophysically relevant surfaces of silicates, carbon 
and ice has been studied by laboratory experiments during the past decade 
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47|, 



51] . From the results of these experiments one can extract, for each surface, 
the energy barriers defined above and the density of adsorption sites. In the simulations 
presented below, the density of adsorption sites was s = 5x 10 13 (sites cm -2 ), the activation 
energies were Eq = 22 meV for diffusion and E\ = 32 meV for desorption. These parameters 



are rounded values of the experimental results for hydrogen recombination on silicates [52 ]. 
The grain temperature was taken as T = 10K. Here we assume that diffusion occurs only 
by thermal hopping, in agreement with the experimental results 52j, |53J]. For small grains, 
it is convenient to replace the hopping rate a (hops s _1 ) by the sweeping rate A = a/S, 
which is approximately the inverse of the time it takes for an H atom to visit nearly all the 
adsorption sites on the grain surface (for a more precise evaluation of the sweeping rate see 



Refs. 
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551] ). The rate equation for this reaction takes the form 

rll N\ 

F — W(N) — 2A(N) , (1) 



dt 

where (N) is the average population size of H atoms residing on the grain. The first term 
on the right hand side of Eq. (TjQ) describes the incoming flux of H atoms. The second 
term accounts for the desorption of H atoms from the surface, which is proportional to the 
number of adsorbed atoms. The third term accounts for the depletion in the population of 
adsorbed H atoms due to the recombination process. 

The production rate of H2 molecules on a single grain, is given by 

R = A(N) 2 . (2) 

For simplicity we assume that these molecules desorb from the surface upon formation. 
The rate equation (TjQ) can be solved either analytically or by numerical integration. The 
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steady state solution can be obtained by setting the left hand side of Eq. ([I]) to zero. Under 
steady state conditions, the average population size of H atoms on a grain is 



(N) 



1 /W 
4 \~A 



(3) 



and the formation rate of molecules per grain is 



R 



ISA 



(4) 



-1 + ^/1 + 87 

where 7 = FA/W 2 . Note that 7 is independent of the grain size. Under these conditions, 
one can define the recombination efficiency 



V 



R 

T/2'' 



(5) 



which is the fraction of adsorbed atoms that come out in molecular form. Note that < 
77 < 1. 

One can identify two limits. In the limit where 7 1 the desorption process is dominant 
and the recombination efficiency is low. In this case, the recombination rate per grain 
is R ~ AF 2 /W 2 and the recombination efficiency is rj ~ 2AF/W 2 . Since R depends 
quadratically on the flux, this regime is referred to as second order kinetics. In the opposite 
limit, where 7 ^> 1, the recombination process is dominant and the efficiency is rj ~ 1. This 
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56|. 



is the regime of first order kinetics, in which R depends linearly on F 

Under given flux and surface temperature, for grains that are large enough to hold many 
H atoms, Eq. (pQ) provides a good description of the recombination process. However, in 
the limit of small grains and low flux, (N) may be reduced to order unity or less. Under 
these conditions, Eq. (JTJ becomes unsuitable, because it ignores the discrete nature of the 
population of adsorbed atoms and its fluctuations. 



B. The Master Equation 

To account correctly for the recombination rate on small grains, simulations using the 
master equation are required. The dynamical variables of the master equation are the 
probabilities P{N) of having a population of N hydrogen atoms on the grain. In the case 
of hydrogen recombination, the master equation takes the form 
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P(N) = F [P(N - 1) - P(N)} + W[(N + 1)P(N + 1) - NP(N)} 

+ A[(N + 2)(N + 1)P(N + 2) - N(N - 1)P(N)], (6) 

where N = 0,1, 2,.... The first term on the right hand side of Eq. (jSJ) describes the effect 
of the incoming flux. The probability P(N) increases when an H atom is adsorbed by a 
grain that already has N — 1 adsorbed H atoms, and decreases when it is adsorbed on a 
grain with N atoms. The second term accounts for the desorption process. The third term 
describes the recombination process. The recombination rate is proportional to the number 
of pairs of H atoms on the grain, namely iV(iV — l)/2. Therefore, the H 2 production rate 
can be expressed in terms of the moments of P{N) as 

R = A((N 2 )-(N)). (7) 

In numerical simulations the master equation must be truncated in order to keep the number 
of equations finite. A convenient way to achieve this is to assign an upper cutoff jV max on 
the population size. The number of equations is thus jV max + 1. The truncated master 
equation is valid if the probability to have a larger population than the assigned cutoff is 
vanishingly small. Therefore, the upper cutoff should be chosen according to the parameters 
of the simulation. 



C. The Moment Equations 

As noted above, the population size of the adsorbed H atoms is given by the first moment 
of P(N), while the reaction rate can be expressed in terms of the difference between the 
second moment and the first moment. Therefore, a closed set of equations for the time 
derivatives of these first and second moments could directly provide all the information that 
we need about the population size and reaction rates [571 ] . Such equations can be obtained 
from the master equation using the identity 

^4r= E N"P{N) t (8) 
at N=0 

where k is an integer. Inserting P{N) according to Eq. (JS)), one obtains the moment 
equations 
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d ^ -. F- W(N) - 2A((N 2 ) - (iV)) 



dt 
d(N 2 ) 

' 1 - F+(2F + W -4A)(N) + (8A-2W)(N 2 ) -4A(N 3 ). (9) 



dt 

This is a set of coupled differential equations, which are linear in the moments (N k ). Al- 
though we have written the equations only for the first two moments, the right hand side 
includes the third moment for which we have no equation. In order to close the set of equa- 
tions one must express the third moment in terms of the first two moments. Different such 
expressions have been proposed. For example, in the context of birth-death processes the 
relation (A^ 3 ) = (N 2 )(N) was used [58!]. This choice makes the moment equations nonlinear, 
with possible effects on their stability. Another common choice is to assume that the third 
central moment is zero (which is exact for symmetric distributions) and use this relation in 
order to express the third moment in terms of the first and second moments 

Here we set up the closure condition by imposing a highly restrictive cutoff on the master 
equation. The cutoff is set at A^ max = 2. This is the minimal cutoff that still enables the 
recombination process to take place. Under this cutoff, one obtains the following relation 



59j. 



between the first three moments 57] 



(N 3 ) = 3(N 2 ) ~2(N). (10) 
Using this result, one can bring the moment equations [Eq. (Q] into a closed form: 



dt 
d(N 2 ) 



F - W(N) - 2A((N 2 ) - (AO) 

F + (2F + W + 4A)(N) - (2W + 4A)(N 2 ). (11) 



dt 

Numerical integration of these equations provides the first two moments, from which the 
population size and reaction rate are obtained. The steady state solution of the moment 
equations can be obtained by setting the derivatives on the left hand side of Eq. ffTTT) to 
zero. One obtains 

(N)= F(A + W) 
v 1 2AF + WA + W 2 

= FJF + A + W) 
x 1 2AF + WA + W 2 K J 
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Using Eq. ([7]) we find that the recombination rate is 



P 

R= — 
W 



s 2 



1 -I i I S_ 

1 T a/VK T W// 



(13) 



In this system one can identify two characteristic spatial scales, namely Si = a/W and 
5*2 = W/ f 57|]. The limit of small grains is obtained when S is smaller than both Si and 
£2. In this limit the populations size of atoms on the grain is small, the rate equation fails 
while the moment equations are accurate. The limit of large grains is obtained when S is 
larger than both Si and S%. In this limit the population size of atoms on the grain is large, 
the rate equation applies and the reaction rate satisfies R oc S. Since F = fS and A = a/S, 
the limit of large grains can be expressed by 



F > W > A. (14) 

We recall that the moment equations were derived by imposing a strict cutoff on the 
master equation. One may thus expect that the equations will apply only in the case of 
small grains, where this cutoff is suitable. However, it turns out that the moment equations 
are valid much beyond this limit. Below we show that the reaction rate obtained from the 
moment equations is accurate for both small and large grains. 

In Fig. [1] we present the population size of H atoms and the reaction rate, vs. grain 
diameter, obtained by the moment equations (circles). The flux of H atoms was F = 
1 x lO -11 ^ (atoms s _1 ). The parameters used in Fig. [1] satisfy the conditions for second- 
order kinetics, namely 7 -C 1. The results are in perfect agreement with the master equation 
(solid lines) for the entire range of grain sizes. For small grains, the rate equation (dashed 
lines) over-estimates the recombination rate, but coincides with the master equation and 
the moment equations for large grains. In Fig. [2] we present the population size of H atoms 
and the reaction rate, vs. grain diameter, under the conditions of first-order kinetics, where 
7 ^> 1. These conditions were obtained by increasing the flux to F = 2.75 x 10~ 8 S (atoms 
s _1 ). The moment equations still provide accurate results for the reaction rate on grains of 
all sizes. However, for large grains the population size of adsorbed H atoms, obtained from 
the moment equations, deviates from the master equation results. 

In conclusion, one can identify two kinetic regimes (first and second order) and two limits 
(small and large grains). In the limit of small grains, the moment equations are valid, as 
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expected, for both kinetic regimes (Table I). In the case of second order kinetics, the moment 
equations are valid both for small grains and for large grains. They provide accurate results 
both for the average number of atoms on a grain and for for the reaction rate. In the case of 
large grains under conditions of first order kinetics the situation is different. The moment 
equations still provide accurate results for the reaction rate but are not suitable for the 
evaluation of the population size on a grain. In the next section we present a more complete 
analysis of the moment equations and their validity. 

D. The Validity of the Moment Equations 

The moment equations were derived on the basis of a very strict cutoff, allowing at most 
two atoms to simultaneously reside on the surface of the grain. Nevertheless, it turns out 
that the equations apply even under conditions in which the population size of adsorbed 
atoms is large. Here we examine the validity of the moment equations in the regimes of first 
and second order kinetics and in the limits of small and large grains. For sufficiently small 
grains the population of adsorbed atoms on a grain is small and is consistent with the strict 
cutoff imposed in the derivation of the moment equations. Therefore in the limit of small 
grains the approximation of Eq. (flUj) is justified, and the moment equations are valid. 

In the limit of large grains the population of adsorbed atoms on a grain is large and 
the rate equation becomes accurate. Therefore, we will test the validity of the moment 
equations in this limit by comparison with the rate equations. In the limit of large grains 
the relations F 3> W and W 3> A are satisfied. We will first consider the case of second 
order kinetics, where FA <C W 2 . Using the relations above, we find that Eqs. (fl2l) are 
reduced to (N) ~ F/W and (iV 2 ) — (N) ~ (F/W) 2 . Thus the results of the moment 
equations for the population size and for the reaction rate coincide with the rate equations. 

In the case of first order kinetics FA 3> W 2 . In the limit of large grains, we obtain from 
Eq. (j!2p that (N 2 ) — (N) ~ Fj (2A), which is consistent with the H 2 production rate obtained 
from rate equations. On the other hand, the expression for the first moment is reduced to 
(N) ~ W/(2A), which is inconsistent with the rate equations, in which (N) ~ yj ' F/(2A). 
These results are in line with the numerical results shown above, indicating that the moment 
equations are suitable for the evaluation of reaction rates on large grains both in first and 
second order kinetics. However, they provide the correct population size of atoms on large 
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grains only in the case of second order kinetics and not in first order kinetics. 

In Fig. [3](a) we show the first moment, (N), as obtained from the moment equations 
(circles) vs. 7, for a very large grain of diameter d = 1mm, where S ~ 1.6 x 10 12 . The results 
are compared with those of the rate equation (dashed line), which under these conditions 
is accurate. For 7 1 (second order kinetics) the rate equation and the moment equations 
coincide. For 7 ^> 1 (first order kinetics) the moment-equation results for the average 
population size deviate from those of the rate equation. Note that even when 7 -C 1, where 
the moment equations apply, the population size of atoms per grain is very large, much 
beyond the cutoffs used in the construction of the equations. In Fig. Mjo) we present the 
production rate of H 2 molecules per grain vs. 7, as obtained from the moment equations 
(circles). The results are in perfect agreement with those of the rate equations (dashed line) 
in the regime of first order kinetics as well as in the regime of second order kinetics. These 
results can be generalized to more complex networks. The analysis becomes more tedious, 
but the conclusions about the domain of validity of the moment equations remain the same. 
From a physical perspective, in first order kinetics the reaction is dominant and desorption 
is suppressed, while in second order kinetics desorption is dominant and the reaction rate is 
low and diffusion limited. 



III. COMPLEX REACTION NETWORKS 
A. The Water-Producing Network 

To generalize the analysis beyond the case of a single species, consider a simple chemical 



network that involves three reactive species: H and O atoms and OH molecules [39l. |4Q|. |60|. 
For simplicity we denote the reactive species by X\ = H, X 2 = O, X 3 = OH, and the 
resulting non-reactive species by X4 = H2, X5 = O2, X$ = H2O. The reactions that take 
place in this network include H + O -> OH (X x + X 2 -> X 3 ), H + H -> H 2 (Xi +X ± -> X 4 ), 
O + O -> 2 (X 2 + X 2 X 5 ), and H + OH -> H 2 (X x + X 3 -> X 6 ). A graph of that 
network is displayed in Fig. H](a). 

Now, the desorption rates of atomic and molecular species on the grain are given by Wi = 
z/-exp[— E\(i) /ksT], where Ei(i) is the activation energy barrier for desorption of species Xi 
and T (K) is the surface temperature. The hopping rate of adsorbed atoms between adjacent 
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sites on the surface is dj = v ■ exp[— Eq(i) /ksT], where Eq(i) is the activation energy barrier 
for hopping of Xi atoms (or molecules). The sweeping rate will be A{ = ai/S. 

The master equation provides the time derivatives of the probabilities P(Ni, N 2 ,N 3 ) that 
a randomly chosen grain will have iVj atoms or molecules of the reactive species X^. It takes 
the form 



P{N U N 2 , N 3 ) = £ F< [P(~, Ni - 1, ..) - P(N U N 2 , N 3 )} 

3 

+ £ W i K N i + Ni + 1,..)- NiP(Ni, N 2 , N 3 )} 

i=i 

+ £ M(Ni + 2)(N { + 1)P(., Ni + 2,..)- Ni(Ni - l)P(N u N 2 , N 3 )} 
i=i 

+(Ai + A 2 ) [{N x + 1)(N 2 + l)P(iV 1 + 1, N 2 + 1, iV 3 - 1) - ^N^Nt, N 2 , N 3 )} 
+(A X + A 3 ) [(JVx + l)(iV 3 + l)P(iV 1 + 1, jV 2 , iV 3 + 1) - iViJV 3 P(JVi, iV 2 , iV 3 )]. (15) 

The terms in the first sum describe the incoming flux, where F{ (atoms s _1 ) is the flux 
'per grain of the species Xi. The second sum describes the effect of desorption. The third 
sum describes the effect of diffusion-mediated reactions between two identical atoms and 
the last two terms account for reactions between different species. The rate of each reaction 
is proportional to the number of pairs of atoms of the two species involved, and to the 
sum of their sweeping rates. The average population size of the X, species on the grain is 
(Ni) = E Ni ,n 2< n 3 N P(Nx, N 2} N 3 ), where N { = 0, 1, 2 . . ., and % = 1, 2 or 3. The production 
rate per grain R(Xk) (molecules s _1 ) of Xk molecules produced by the reaction Xi+Xj — > X^ 
is given by R(X k ) = (A { + A j ){N i Nj), or by R(X k ) = A^N^Ni - 1)) in case that % = j. 

For complex networks with more than one species the truncation of the master equation 
demands setting upper cutoffs for all the reactive species, iV™ ax , i — 1, . . . , J, where J is the 
number of reactive species. The number of coupled equations is then Ne = n/=i(^ r l max + 1), 
which grows exponentially with the number of reactive species. This severely limits the 



applicability of the master equation to interstellar chemistry [43l. |44|] . To reduce the number 
of equations one tries to use the lowest possible cutoffs under the given conditions. In any 
case, to enable all reaction processes to take place, the cutoffs must satisfy A^ max > 2 for 
species that form homonuclear diatomic molecules (H2, O2, etc.) and N^ ax > 1 for other 
species. 
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As in the case of hydrogen recombination, the population sizes and reaction rates are 
completely determined by all the first moments and selected second moments of the distri- 
bution P(Ni, N 2) N 3 ). Therefore, a closed set of equations for the time derivatives of these 
first and second moments could provide the complete information on the population sizes 
and reaction rates. For the simple network considered here one needs equations for the time 
derivatives of the first moments (iVj), i — 1, ... ,3 and of the second moments (Nf), (iVf), 
(NiN 2 ) and (N1N3). We obtain these equations from the master equation using the identity 



TCI Affe ATC\ J \ 



d{N;N 2 Ni) = ^ Na^cp^N^ (16) 

0,1 Ni,N 2 ,N 3 =0 

where a, b, c are integers. Here we show three of the resulting moment equations: 



^ Ai) = Fl - Wl (N l )-2A 1 ((N^)-(N 1 ))-(A 1 + A 2 )(N l N 2 )-(A 1 + A 3 )(N l N 3 ) 

f// = F 1 + (2F X + Wi- AAJiNi) - (2W 1 - 8A 1 )(N^} - AA 1 (N^) 

+ {A 1 + A 2 )((iV 1 iV 2 ) - 2(^2)) + {A 1 + A 3 )((AyV 3 } - 2(N*N 3 )) 

F!{N 3 ) + F 3 (N ± ) - {W l +W 3 - 3A, - A 3 )(NiN 3 ) - (3A, + A 3 )(N*N 3 ) 



dt 
d(Nf) 



dt 

- (A, + A 3 )(iViiV3 2 ) " ( A i + ^ 2 )((iViiV 2 ) - (NfN 2 ) + (N.N.N,)). (17) 

As expected, the right hand sides of these equations include third order moments for which 
we have no equations. If we add equations for their time derivatives, they will include fourth 
order moments, and again will not enable us to close the set of equations. In order to close 
the set of moment equations we must express the third order moments in terms of the first 
and second order moments. We recall that in the case of a single reactive specie this was 
achieved by imposing an extremely restrictive cutoff on the master equation. Generalizing 
it to our present system we choose the minimal cutoffs that do not terminate any of the 
reactions taking part in the network. These cutoffs allow at most two atoms or molecules 
to be adsorbed on the surface at any given time. Furthermore, any two atoms or molecules 
that reside simultaneously on the surface must belong to species that react with each other. 
In our network these cutoffs allow only eight non- vanishing probabilities, namely, P(0, 0, 0), 
P(0, 0, 1), P(0, 1,0), P(l, 0, 0), P(2, 0, 0), P(0, 2, 0), P(l, 1, 0) and P(l, 0, 1). Expressing the 
third moments, with these cutoffs, leads to the following results: 
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(N?) = 3(JV?>-2(JVi> 
(N?Nj) = {NiNj) 
(NiNjNk) = 0. (18) 



Using these rules to modify Eqs. ffTTj) one obtains a closed set of the form 
d(Nx) 



dt 
d(N 2 ) 

dt 
d(N 3 ) 

dt 
d(N?) 
dt 

d(N 2 2 ) 
dt 

dt 

djNiNs) 
dt 



Fx - Wx(Nx) - 2A 1 ((Nf) - (Nx)) - (Ax + A 2 )(NxN 2 ) - (Ax + A 3 )(NxN 3 ) 

F 2 - W 2 (N 2 ) - 2A 2 ((Nl) - (N 2 )) - (Ax + A 2 )(NxN 2 ) 

F 3 - W 3 (N 3 ) - (Ax + A 3 )(NxN 3 ) + (Ax + A 2 )(NxN 2 ) 

Fx + (2Fx + Wx+ 4Ax)(Nx) - (2Wx + 4A 1 )(Ar i 2 ) 
(Ax + A 2 )(NxN 2 ) - (Ax + A 3 )(NxN 3 ) 

F 2 + (2F 2 + W 2 + AA 2 )(N 2 ) - (2W 2 + ±A 2 )(N 2 2 ) - (Ax + A 2 )(NxN 2 ) 
Fx(N 2 ) + F 2 (Nx) - (Wx + W 2 + Ax + A 2 )(NxN 2 ) 

Fx(N 3 ) + F 3 (Nx) - (Wx + W 3 + Ax + ^) (JV X JV 3 ) . (19) 



This set of equations is the minimal set required for the system. It includes one equation 
that accounts for the population size of each reactive specie and one equation that accounts 
for the reaction rate of each reaction. As in the hydrogen system, the equations were derived 
using extremely strict cutoffs, that are expected to apply only in the limit of small grains and 
low flux. Nevertheless, the moment equations are valid well beyond the restrictive cutoffs 
under which they were derived. In Fig. [5]^a) we present the population sizes of H ( x ) and O 
(+) atoms on a grain vs. grain diameter, obtained from the moment equations. In Fig. GO^b) 
we show the production rates of H 2 (circles), 2 (triangles) and H 2 (squares), obtained 
from the moment equations. The results are in excellent agreement with the master equation 
and in the limit of large grains they also coincide with the rate equations. The parameters 
used in the simulations are Fx = 5.0 x lO -10 ^ (atoms s -1 ), F 2 = O.LFi and F 3 = 0. The 
activation energies for diffusion and desorption were taken as E (l) = 44, Ex(l) = 52, 
E (2) = 47, Ex(2) = 54, E (3) = 47 and ^(3) = 54 meV for H, O and OH respectively. 
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The surface temperature of the grain was T = 15K. For the species other than H there are 
no concrete experimental results, and the values reflect the tendency of heavier species to 
bind more strongly. 

The generalization to more complex networks is straightforward. As explained above, 
the cutoffs imposed on the master equation give rise to a minimal set of non-vanishing 
probabilities, that allow all the reactions to take place. These non-vanishing probabilities 
are P(0, . . . , 0, iV» = 1, 0, . . . , 0) for the each node, i, and P(0, . . . , 0, N t = 1, 0, . . . , 0, N, = 
1, 0, . . . , 0) for each edge connecting nodes i and j. A loop connected to the node i corre- 
sponds to the probability P(0, . . . , 0, Ni = 2, 0, . . . , 0). All other probabilities vanish. In the 
realm of two-body reactions additional probabilities such as P(3, 0, . . . , 0) or P(2, 1, 0, . . . , 0) 
are never required. Under these restrictions, one can easily show that all the third moments 
that appear in the moment equations can be expressed in terms of second moments according 
to Eq. (HB|). 



B. The Methanol Network 



Let us now examine how the moment equation method performs in a more complex 
network. Recent laboratory experiments provide evidence that methanol in interstellar 



clouds is formed via reaction networks on dust grains |35l. |36| . To simulate these networks, 
consider the case in which a flux of CO molecules is added to the network shown in Fig. H](a). 
This gives rise to the network shown in Fig. H^b), which includes the following sequence of 



hydrogen addition reactions 



43|: H + CO -> HCO (X 1 +X 7 ^ X 8 ), H + HCO -> H 2 CO 



(X! + X 8 -> X 9 ), H + H 2 CO -> H3CO {X 1 + X 9 -> X10) and H + H 3 CO -> CH 3 OH 
(Xi + X 10 — > An). Two other reactions that involve oxygen atoms also take place: O + 
CO -> C0 2 (A 2 + A 7 -> A 12 ) and O + HCO — > C0 2 + H (A 2 + A 8 -> A 12 + X x ). This 
network was studied before using the multiplane method, which required about a thousand 
equations compared to about a million equations in the master equation with similar cutoffs. 
The moment equations include one equation for each species (node in the graph) and one 
equation for each reaction (edge in the graph), namely, the network shown in Fig. H^b) 
requires only 17 equations. We have performed extensive simulations of this network using 
the moment equations and found that they are in agreement with the master equation results. 
In Fig. Efa) we present the population sizes of H (+), O (*) and CO (x) on a grain vs. 
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grain diameter, as obtained from the moment equations. In Fig. Mjo) we present the results 
of the moment equations for the formation rates of CH3OH (squares), H2O (circles) and 
O2 (triangles) and compare them to the master equation (solid line) and the rate equations 
(dashed lines). As before, the moment equations prove to be consistent with the master 
equation, and for large grains with the rate equations as well. In the simulations presented 
here the activation energies for diffusion and desorption of CO, HCO, H2CO and H3CO 
were taken as E (7) = 50, E 1 {1) = 55, E (S) = 52, £1(8) = 58, E (9) = 53, £i(9) = 59 
and .E'o(lO) = 55, E%(10) = 62 meV, respectively. The flux of CO molecules was taken as 
F 7 = 0.2Fi. 

In complex reaction networks the reduction in the number of equations becomes extremely 
significant. For each reactive species, X, (denoted by a node in the graph), there is one 
equation for the time derivative of the first moment (N{). For each reaction, between species 
Xi and Xj (represented by an edge), there is one equation for the time derivative of the 
second moment (NiNj). The number of equations is thus j\f^ oment — N n + N c , where N n and 
N e are the numbers of nodes and edges, respectively. Most reaction networks are sparse, 
namely the number of edges is of the same order of magnitude as the number of nodes. As 
a result, the number of moment equations is roughly linearly dependent on the number of 
species. This is in contrast with the exponential growth of the master equation and the 
polynomial growth of the multiplane equations jlt^ . 

IV. DIAGRAMMATIC CONSTRUCTION OF THE EQUATIONS 

In the previous Section we described a procedure for the construction of the moment 
equations, that consists of three steps. In the first step the desired set of moment equations 
is determined according to the network graph. This set includes one equation for the first 
moment associated with each species (node) and one equation for the second moment as- 
sociated with each reaction (edge or loop). In the second step the time derivative of each 
of these moments is expressed as a suitable linear combination of first, second and third 
moments. In the third step the equations are brought into a closed form, by expressing all 
third moments in terms of first and second moments according to Eq. (JTBl) . This procedure 
is suitable for simple networks but becomes difficult to execute as the network complexity 
increases. 
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Here we present a simple and highly efficient diagrammatic approach for the construction 
of the moment equations for any given network. Consider the reaction networks shown in Fig. 
HI Each node in the graphs represents one of the reactive species, and each edge stands for 
a reaction between two species. A loop represents a reaction between two atoms/molecules 
of the same species. The list of moments that should be included in the moment equations 
is easily obtained from the graph: (a) For each each node, i, there is an equation for the 
first moment (N); (b) For each edge connecting the nodes i and j there is an equation for 
the second moment (NNj); (c) for each loop attached to node % there is an equation for the 
second moment (Nf) . 

To derive the moment equations from the graph of a given network, we show how each 
element in the graph is expressed in the master equation and then translate it into a cor- 
responding term in the moment equations. Consider first the contribution of nodes. Each 
node represents a species, that may be adsorbed on the surface or desorbed from it. In the 
master equation, these two processes are described by 

Q _^ F i [P(...,N i -l,...)-P(...,N i ,...)} (20) 
+ WiUNi + 1)P(. . . , Ni + 1, . . .) - NP(. ..,Ni,...)\. 

To construct the equation for (Ni) we sum over all the terms of the master equation, ac- 
cording to Eq. f|T6|) . The terms above then make the following contribution to the moment 
equations 

© -^Fi-WiiNt). (21) 

In case that the species i reacts with itself, one needs to include an equation for (Nf). As 
before, a proper summation over the master equation is performed, but here the contribution 
of the node % will be 

© — ► Fi + 2Fi(N) - 2W l (Nf) + Wi(Ni). (22) 

If the species i reacts with some other species, j, it will be necessary to include an equation 
for (NiNj). The contribution of the node i to this equation is 

© ^F^-W^NNj). (23) 
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Note that in the equation for (NNj), one should also include the analogous term for the 
node j. 

Consider the contribution of edges to the moment equations. Each edge represents a 
reaction between two different species. In the master equation, the term related to the edge 
between i and j is 

(Ai + A 3 ) [(iVj + 1)(A^ + 1)P(. ..^ + 1,^ + 1...)- A^P(. . . N t , Nj ...)]. 

(24) 

The contributions of this term to the moment equations appear in all the equations for 
moments involving either species i or species j. For example, the contribution to the equation 
for (NiNj) is 

© 

Y — -(A i + A j )((NfN j ) + (N i Nf)-(N i N j )). 

© (25) 

After applying the rules of Eq. (1TB]) . this term simplifies to — (Ai + Aj)(NiNj). 

In general, each moment equation is obtained by a proper summation over terms of the 
master equation, using Eq. (|T6l) . Each term in the master equation describes a certain 
process such as adsorption, desorption or reaction. Consider the equations for (iVj) and 
(Nf). When the summations are carried out, the only non- vanishing contributions are from 
processes that involve the species Xj. Similarly, consider the equation for (NiNj). In this 
case the only non-vanishing contributions are from processes that involve either the species 
Xi or the species Xj. This gives rise to a dramatic simplification in the construction of the 
moment equations. In the equation for any desired moment we need to consider only those 
master equation terms related to the species that appear in that particular moment. In 
the diagrams, this means including only graph elements that are directly connected to the 
corresponding nodes. The complete set of translation rules from a given graph element to 
the corresponding terms in the moment equations is given in Tables UTI - IIVI 

To demonstrate the construction of the moment equations we reconsider the reaction 
network shown in Fig. H^a). First, we identify the relevant moments associated with the 
nodes, edges and loops in the network graph. In this particular network, the nodes call for 
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(iVi), (N 2 ) and (N 3 ). The edges call for (NiN 2 ) and (iViiV 3 ). Finally, the two loops in the 
graph require the inclusion of (Nf) and (iV|). Altogether, there are seven equations. Let us 
begin with the equation for (Aq). This species is represented by a node connected to two 
edges and one loop. Using our diagrammatic scheme, this equation takes the form 




This equation simply consists of all the graph elements related to the species X\. That 
diagrammatic sum can be translated, using Table [Til giving rise to the equation 

(Aq) = iq - W 1 (N 1 ) - (A 1 + A 2 )(N 1 N 2 ) 

- (A, + A 3 )(N 1 N 3 ) - 2A 1 ((A 1 2 > - (Aq)) (27) 

The same diagrammatic equation [Eq. (126]) ] is used when writing the equation for the time 
derivative of (Nf). However, the translation rules from graph elements to equation terms 
are different and are given in Table IIHI Using these rules we obtain 

(N?) = F x + (2F 1 + W 1 + 4A 1 )(A 1 ) - {2W t + 4A 1 )(Ar i 2 ) 

- (A 1 +A 2 )(N 1 N 2 )-(A 1 +A 3 )(N 1 N 3 ). (28) 

Next, we construct the moment equation for (NiN 2 ). This time we need to include some 
additional graph elements, namely, those related to X 2 , so the diagrammatic equation takes 
the form 

(NiN 2 ) = © + (D + T+T+0 + Q- ( 29 ) 

© (3) ® ® 

This is an equation for a mixed second moment, that should be translated according to 
Table [TV] The resulting equation is 

(AyV 2 ) = Fi(N 2 ) + F 2 (N 1 ) - (W 1 + W 2 ){N l N 2 ) - (A 1 + A 2 )(N 1 N 2 ). (30) 
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The last equation that we demonstrate explicitly is for the moment (N 3 ). This time the 
diagrammatic equation will include an additional graph element stating the fact that X 3 is 
produced by a reaction between the other two species in the network. The diagrammatic 
equation is 




which, according to Table [TT] translates to 

(N 3 ) = F 3 - W 3 (N 3 ) - (A, + A 3 )(i\qAr 3 ) + (A + A 2 )(N 1 N 2 ). (32) 

This procedure enables an efficient and straightforward construction of the set of moment 
equations for any given graph, regardless of its complexity. The construction process can be 
easily automated into a computer program that receives the network structure as an input 
and provides the set of equations as an output. This opens the way to efficient stochastic 
modeling of complex reaction networks of any size. 

V. SUMMARY 

We have presented an efficient method for the stochastic simulation of complex reaction 
networks. The method is based on moment equations. It accounts correctly for the reaction 
rates not only for macroscopic populations, where rate equations apply, but also in the 
limit of low copy numbers where fluctuations dominate. The method exhibits several crucial 
advantages over existing techniques for stochastic simulations. The number of equations 
is reduced to the absolute minimum in a stochastic simulation, namely one equation for 
each reactive species and one equation for each reaction. For any given network the set of 
moment equations can be constructed easily and efficiently using a diagrammatic approach. 
The equations are linear, thus, for steady-state conditions they can be easily solved using 
algebraic methods. 

We have demonstrated the method for reactions taking place on dust-grain surfaces in 
interstellar clouds. We have shown that it applies very well even under the extreme inter- 
stellar conditions of low gas density and sub-micron grain sizes. The moment equations can 
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be easily coupled to the rate equation models of interstellar gas-phase chemistry. Therefore, 
the proposed method is expected to enable the incorporation of grain-surface chemistry into 
these models. 

We expect this method to be useful for the simulation of many other systems that exhibit 
a related mathematical structure [61]. One such application is for reactions takin g p lace 



on terrestrial surfaces such as metal catalysts in nanometric systems 
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231 ] . Another example is aerosol chemistry in stratospheric clouds, where complex reaction 



networks involving aerosol particles take place 62|. The gas density is not nearly as low 
as in the interstellar clouds. However, the small size of the aerosol particles and the low 
density of some of the reactive species may require stochastic methods. Another example is 



genetic networks in cells and bacteria 
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64], 



65. 



66j . These networks describe the process of 



protein synthesis and its regulation. They include the interactions between genes, mRNA's 
and proteins. Since some of these components appear in low copy numbers, the simulation 
of these networks requires stochastic methods. 

We thank Azi Lipshtat for helpful discussions. This work was supported by the Israel 
Science Foundation and the Adler Foundation for Space Research. 
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TABLE I: The moment equations method is valid for producing the first and second moments of 
the distribution P(N) both for large and small grains. These moments enable the evaluation of 
the population size of atoms per grain, and the production rate of H2 molecules on a grain. The 
one exception is for producing the first moment for large grains in the regime of first order kinetics. 
In that case, one may simply use the rate equation model, which is valid for large grains, or the 
corrected moment equations which always apply for large grains. 







The Validity of the Moment Equations 
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TABLE II: Diagrammatic approach for the construction of the moment equations. The equation 
for (Ni) is given by the sum of all graph elements related to the species X{. To construct this 
moment equation, each graph element is substituted by the Reduced Equation Term, obtained 
using the rules of Eq. (I18p . The middle column provides the exact terms of the moment equation, 
before the approximation of Eq. (fT8|) is applied. 



Graph Elements Translation for the (Ni) Equation 
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Ak((N 2 k ) - (N k )) 


A k ((Nl) - (N k )) 



26 



TABLE III: The equation for (Af) is the sum of all graph elements related to the species Xj. 
Some of the exact equation terms, given in the middle column, include third order moments. The 
reduced equation terms, given in the third column, consist of only first and second order moments. 



Graph Elements Translation for the (N?) Equation 
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TABLE IV: The equation for (NiNj) is the sum of all graph elements related to the species 
and Xj. Elements related to both species should be included only once. 
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Graph Elements Translation for the (NiNj) Equation 
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FIG. 1: The population size of H atoms (a) and the production rate of H2 molecules (b) on the 
surface of a grain vs. grain diameter, d, under conditions of second-order kinetics. The number 
of adsorption sites, S, is also shown. The results of the moment equations (circles) are in perfect 
agreement with the master equation (solid line) for all grain sizes. In the limit of large grains they 
also coincide with the rate equation (dashed line). 
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FIG. 2: The population size of H atoms (a) and the production rate of H2 molecules (b) on the 
surface of a grain vs. grain diameter, d, under conditions of first-order kinetics. The production 
rate obtained by the moment equations (circles) is in perfect agreement with the master equation 
(solid line) for all grain sizes. However, the population size deviates for large grains. The rate 
equation results are also shown (dashed lines). 
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FIG. 3: (a) The population size of H atoms per grain vs. the parameter 7, as obtained from the 
moment equations for a very large grain of diameter d = 1mm (circles). Under these conditions 
the rate equation is accurate. For 7 < 1 (second order kinetics), the moment equations are in 
perfect agreement with the rate equation (dashed line). This is in spite of the large population 
size. For 7 > 1 (first order kinetics), the moment equations do not provide accurate results for the 
population size; (b) The production rate of H2 molecules per grain vs. 7, as obtained from the 
moment equations (circles). The results are in perfect agreement with those of the rate equation 
(dashed line) both for second order and for first order kinetics (even when there is a deviation in 
the first moment). 
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FIG. 4: (a) The reaction network on a grain exposed to fluxes of H and O atoms and OH molecules. 
Each node represents a reactive species and each edge represents a reaction. The reaction products 
are specified near the edges; (b) The network obtained by adding CO molecules, which gives rise 
to the production of methanol. 
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FIG. 5: (a) The population sizes of H (x) and O (+) atoms on a grain vs. grain diameter, d, 
obtained from the moment equations for the reaction network presented in Fig. 3(a). The results 
are in excellent agreement with the master equation (solid line). The rate equation results are also 
shown (dashed line), (b) The production rates of H2 (circles), H2O (squares) and O2 (triangles) 
per grain vs. d, obtained by the moment equations. The results are in excellent agreement with the 
master equation (solid line) . In the limit of large grains they also coincide with the rate equations 
(dashed line). Note that the moment equations are valid both in the limit of small population sizes 
(when d is small), and even when the population is high (large ds). The number of adsorption 
sites, S, is also displayed. 
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FIG. 6: (a) The population sizes of H (+), O (*) atoms and CO molecules (x) on a grain vs. grain 
diameter, d, obtained from the moment equations for the reaction network presented in Fig. 3(b). 
The results are in excellent agreement with the master equation (solid line). The rate equation 
results are also shown (dashed line), (b) The production rates of H2O (circles), CH3OH (squares) 
and O2 (triangles) molecules on a grain vs. d, obtained from the moment equations. The results 
are in agreement with the master equation (solid line). For small grains, the rate equation results 
(dashed line) show significant deviations. 
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